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This  work  establishes  the  feasibility  of  handling  problems  of  projectile/bore 
interaction  via  the  method  of  finite  elements.  The  problem  treated  herein 
involves  the  analytical  determination  of  the  time-variant  matrix  coefficients 
for  the  structure  equations  of  motion  in  a manner  consistent  with  the  finite 
element  displacement  method.  The  resultant  equations  of  motion  are  solved 
numerically  and  results  compare  favorably  with  those  reported  by  experi- 
mentalists in  the  literature.  For  this  purpose,  the  otherwise  general 
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procedure  is  applied  to  the  problem  of  predicting  the  response  of  a simply 
supported  beam  to  a concentrated  mass  in  motion  at  constant  velocity,  Tli* 
technique  can  be  easily  extended  to  gun  tube  geometries  by  merely  altering 
the  boundary  conditions,  increasing  the  number  of  finite  elements  and  allowing 
their  dimensions  to  approximate  the  tube  of  interest.  Similarly,  the  feature 
of  variable  projectile  velocity  can  easily  be  included. 
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SUMMARY 


The  results  displayed  in  this  report  demonstrate  the  feasibility  of 
handling  moving  mass  problems  by  the  finite  element  displacement  method. 
There  is  no  apparent  reason  why  the  method  could  not  be  applied  to  gun 
tube  geometries  - should  engineering  fund  allocations  be  made  available 
for  this  purpose.  (Nearly  all  of  the  work  contained  in  this  report  was 
funded  through  basic  in-house  research  allocations.) 

In  the  interest  of  inproved  efficiency,  more  detailed  study  should 
be  given  to  the  choice  and  use  of  integration  algorithms.  In  the  IBM 
version  of  Hamming's  method  (employed  in  this  work),  the  user  is  allowed 
certain  freedom  to  choose  error  weights.  This  choice  can  affect  ac- 
curacy and  run  times  and  therefore  should  be  arrived  at  on  a rational 
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INTRODUCTION 


Within  the  U.S.  Army  there  has  always  been  an  interest  in  the  bal- 


listic response  of  the  gun  tube  and  projectile  during  firing.  This 


interest  has  several  justifications  - not  the  least  of  which  are  the 


effects  of  tube  motion  on  round  accuracy  and  the  reactions  inuucea 


between  the  interior  bore  surface  and  the  projectile  - the  latter  hav- 


ing to  do  with  the  structural  integrity  of  the  projectile  casing  and 


the  wear  caused  by  excessive  projectile/bore  interface  pressures. 


During  1974  this  laboratory  reported  a NASTRAN  finite  element  anal- 


ysis of  the  axi3ymmetric  response  of  an  unconstrained  M113  gun  tube  to 


the  ballistic  pressure/time  loading  corresponding  to  a 'zone  three' 


charge.  Also  computed  in  detail  was  the  transverse  response  to  the 


curvature-induced  'Bourdon'  load.  Less  detailed  calculations  alluded 


to  the  potential  response  which  could  be  expected  from  the  tube  as  a 


result  of  the  moving  concentrated  weight  of  the  projectile  along  the 


bore  surface.  Finally,  hope  was  expressed  for  future  progress  in  hand- 


ling the  various  aspects  of  the  moving  projectile  mass. 


The  Importance  of  the  Moving  Mass  Problem 


In  more  conventional  applications  the  difference  in  the  response  of 


a structure  to  a moving  mass  and  to  the  weight  of  this  mass  is  not  very 


great.  However,  if  the  moving  mass  is  an  appreciable  fraction  of  the 


total  structural  mass  of  the  problem  - or  if  the  velocity  of  the  mass 


is  very  large,  an  unconventional  structural  problem  is  defined  in  which 


the  difference  may  be  quite  pronounced.  There  ave  at  least  two  illustra 


I\1 


tions  which  come  to  mind:  (1)  - a railroad  train  crossing  a long 

trestle* (in  which  the  moving  mass  compares  to  the  mass  of  the  trestle) 
and  (2)  - a nrojectile  traversing  a long  tube  (in  which  uncommonly 
high  velocities  are  evident).  We  are  professionally  interested  in 
the  latter  as  a knowledge  of  the  response  of  a tube  to  a moving  mass 
leads  immediately  to  a good  estimate  of  the  resultant  interface  pres- 
sure between  the  projectile  and  the  bore  surface  - an  important  first 
step  toward  the  ultimate  understanding  of  the  causes  of  bore  and 
muzzle  wear  problems  and  prooer  design  of  projectile  casings  and 
rotating  bands.  The  importance  of  the  structural  response  of  the 
tube  in  altering  the  round  trajectory  is  also  apparent.  Thus  treat- 
ment of  the  moving  mass  problem  has  a wide  base  of  military  justi- 
fication. 

State  of  the  Art  in  Moving  Mass  Problems 

Although  the  method  of  computation  reported  herein  is  immediately 
applicable  to  a finite  element  beam  model  of  a gun  tube,  it  is  first 
necessary  to  assess  the  accuracy  of  the  mothod  by  applying  it  to  a 
problem  which  has  received  prior  treatment  and  reported  in  the  liter- 
ature. The  most  common  base  for  comparison  has  been  the  experimental 

2 

work  of  Ayre,  Jacobsen,  and  Hsu  - Presumably  because  of  the  lack  of 


(1)  Stokes,  Sir  George  G.,  "Discussion  of  a Differential  Equation  Re- 
lating to  the  Breaking  of  Railway  Bridges,"  Trans  Cambridge  Phil  Soc, 
8,  p.  707,  1849 

(2)  Ayre,  R.S.,  Jacobsen,  L.S.,  and  Hsu,  C.S,,  "Transverse  Vibration 
at' One  and  of  Two  Span  Beami  Under  the  Action  of  a Moving  Mass  Load," 
Prcc.  of  First  National  Congress  on  Applied  Mechanics,  June  1951 
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any  simple  analytical  solution  with  which  to  compare  results,  ( The 
analytical  solution  of  Schallenkamp  is  not  convenient  for  general 
application  - this  will  be  mentioned  further  on  in  the  report) 

This  report  is  therefore  concerned  with  nredicting  the  resnonse 
of  a uniform, simply  sunported  beam  while  subjected  to  a concentrated 
mass  moving  along  its  length  at  constant  velocity  under  the  influence 
of  gravity.  The  method  enmloyed,  however,  is  immediately  applicable 
to  time  variant  mass  velocities  as  well  as  other  boundary  conditions 
and  variable  beam  cross  section. 

Before  commencing  with  the  details  of  the  moving  mass  problem  it 
may  be  worthwhile  to  point  out  the  theoretical  differences  between 
nrohlems  in  v/hich  masses  are  in  motion  and  those  which  involve  only 
moving  forces . 

To  begin  with,  moving  loads*  are  but  special  cases  of  time  and 
space  variant  forcing  functions  f(x,t)  for  one  dimensional  structures 
such  as  beams.  The  customary  beam  equation  of  forced  motion  can  be 
written  3: 


MB  ♦ Iw  ■ f(r.,t) 
where  M and  L are  operators: 

M « oA;  L ■ 


« « « 


0) 


2 2 
(El 3 /ax  ) 


(3)  Tong,'  K„  Theory  of  Mechanical  Vibrations  , J,  Wiley  5 Sons,  I960 
n.  300 


* The  term  'load*  is  intended  to  be  general  and  to  represent  applied 
loads  such  as  moving  forces  and/or  masses 


Solving  the  homogeneous  form  of  equation  (1)  (f(x,t)  » 0),  i.e,,  the 
free  vibration  problem-leads  to  the  eigenfunctions  rn(x)  and  the 
eigenvalues  «n  with  which  any  forced  motion  problem  (f(x,t)  ^ 0) 
can  be  solved. 

Transforming  to  modal  coordinates  p^(t)  defined  in  the  expansion: 

w(x,t)  * Z PjftJr^x)  ...(2) 

leads  to  an  infinite  number  of  differential  equations  which  are 
uncoupled: 

2 

i.e,  »li(Pi  ♦ «iPi)  ■ ^{x,t)ri(x)dx  .,,(3) 

l 

where  ra^  ■ ^rjMCrjJdx 

Now  if  f(x,t)  is  a moving  concentrated  force,  i.e,, 

f(x,t)  - -P06(x-vt) 

then  (3)  becomes: 

®ii (Pi  ♦wjpl)  - -P0rt(vt)  .,.(4) 

The  solution  to  (4)  can  be  written  immediately  in  terms  of  the 
eigenfunctions  and  eigenvalues  of  the  free  vibration  problem. 

For  a moving  mass,  however, 

f(x,t)  ■ -nipW  6 (x-vt) 
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When  substituted  in  (3),  the  right  hand  side  becomes; 


mp^i^w/3t2rjL(x)fi(x-vt)dx  « 32w/3t2ri(x)  | (-m*) 

x»vt 


but  3$/3t2l  * w(vt)  ♦ 2#'(vt)  ♦ v2w"(vt)  . • m 3/3x 

x»vt 


using  (2);  32w/3t2|  ■ £{r,j(vt)pj  ♦ 2vr^(vt)pj  ♦ v2r!J(vt)pj} 


then  (3)  becomes: 

mii(pi  ♦ f (rj(vt)pj  ♦ 2vrj(vt)pj  ♦ v2r»(vt)pj) 

• • • (5) 

In  contrast  to  (4),  equations  (S)  are  not  uncoupled.  In  fact  all  of 

the  variables  pi#  i ■ 1,  • , appear  in  each  of  the  infinite  number 

of  differential  equations.  To  make  matters  worse,  each  variable  has 

a time-variant  coefficient.  To  date  the  only  exact  mathematical 

a 

treatment  appears  to  be  that  due  to  Schallenkamp  involving  a triple 
infinite  series  equation  for  unknown  Fourier  coefficients. 

Thus  the  whole  concept  of  natural  frequencies  and  modes  Of  vi- 
h \tion  loses  its  value  in  quantitative  determination  of  the  response 
of  a system  with  time  variant  properties  • of  which  the  moving  mass 
is  but  a special  case.  For  every  location  of  the  mass  along  the  beam 
we  have  a new  infinity  of  eigenvalue  solutions.  With  an  infinite  num- 

(4)  Schallenkamp,  A.,  "Schwingungen  von  Tragern  bei  bewegten  Lasten," 
Ingenieur-Archiv,  v.8,  1937,  pp.  182-198 
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ber  of  locations  for  the  mass  to  occupy,  we  thus  have  a double  infinity 
of  eigenfunctions  and  eigenvalues. 

Equation  (1)  can,  however,  be  solved  numerically  - regardless  of 
the  space  and  time  dependent  material  properties  and  the  forcing  term. 
In  what  follows  we  eliminate  the  space  variable  through  the  finite 
element  process.  The  time  variable  is  handled  by  a numerical  in- 
tegration procedure  of  common  variety  (predictor-corrector).  In 
essence  therefore,  we  approximate  our  continuum  description  (1),  by 
a finite  number  of  ordinary  differential  equations  with  time  depen- 
dent coefficients. 

Continuum  Description  of  the  Moving  Mass  Problem 

* 

In  previous  work  a derivation  was  tiven  for  the  equation  of 
forced  transverse  motion  of  a beam  model  of  the  M113  gun  tube: 

(Ely")”  ■ -k(x,t)y"  ♦ (pAX0 (t)  ♦ pAgsino)y"(x-l)  ♦ (pAX0(t)  ♦ 

♦ pAgsina)y*  - nip (y  ♦ 2v?’  ♦ gcosa  ♦ v4y")5(x-vt)  - pgAcosa  - pAy 

• • • (6) 

This  equation  assumes  that  A,  the  beam  cross  section  is  uniform  and 

that  v - the  velocity  of  the  moving  mass  mp  is  constant.  In  that  we 

will  be  applying  the  equation  to  beam  elements  of  uniform  and  equal 

cross  section  the  former  assumption  is  consistent.  This  restriction 

is  not  necessary  to  the  generality  of  the  method.  The  assumption  of 

(7)  Simkins,  T.E.,  Pflegl,  G.  and  Scanlon,  R. , "Dynamic  Response  of 
the  Ml 13  Gun  Tube  to  Travelling  Ballistic  Pressure  and  Data  Smoothing 
as  Applied  to  XM15C  Acceleration  Data,"  Watervliet  Arsenal  Technical 
Report  WVT-TR-75015,  April  1975 


✓ ./I 


• •!«••••« 


constant  velocity  is  motivated  by  our  desire  to  compare  results  with 

the  work  of  Ayre,  Jacobsen  and  Hsu,  In  (6)  the  term  k(x,t)y"  is  the 

•• 

so  called  'Bourdon*  force.  Terms  in  X0  represent  transverse  forces 
induced  by  axial  recoil  acceleration  when  beam  curvature  and  slope 
are  non-zero.  Terms  in  pAg  are  due  to  the  beam  weight.  None  of 
these  forces  are  of  interest  in  this  report  which  deals  exclusively 
with  forces  induced  by  the  moving  mass: 

i.e.  f(x,t)  ■ -rap(y  + 2vf'  * g ♦ v2y")6(x-vt)  ,,,(7) 

where  a,  the  tube  elevation  angle,  has  been  made  equal  to  zero. 

Thus  the  special  version  of  equation  (6)  we  will  be  concerned  with 
in  this  report: 

Ely""  ♦ pAjf  ■ -mp(y  ♦ 2vf ' ♦ g ♦ v2y")5(x-vt)  ...(8) 

where  E,I,p,  and  A assume  constant  values.  The  left  hand  side  of 
(8)  is  recognized  as  originating  from  the  simplest  of  beam  theory, 
i.e.,  where  the  entire  transverse  deflection  of  the  beam  is  assumed 
to  be  due  to  bending  moment  only.  The  right  side  of  (8)  therefore 
represents  the  totality  of  applied  loads,  the  first  term  corresponding 
to  the  inertia  of  rr^;  the  second  is  a 'Coriolis'  type  load;  the  third 
is  due  to  the  gravitational  force  on  m^  (its  weight)  and  the  fourth 
is  the  'centrifugal'  force  due  to  nip  following  the  beam  curvature. 

The  Dirac  - delta  function  -8-  specifies  that  each  of  these  forces 
acts  in  a concentrated  fashion  at  the  location  x ■ vt  along  the  beam. 
As  written  equation  (8)  is  a continuum  description  of  the  problem. 
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„ r; 


Using  the  method  of  finite  elements  the  space  variable,  x,  will  be 
discretized  - resulting  in  a set  of  ordinary  differential  equations 
with  time  as  the  independent  variable. 

THE  FINITE  ELEMENT  DISPLACEMENT  METHOD 

The  basic  procedure  in  finite  element  procedures  is  to  consider 
the  structure  of  interest  as  being  composed  of  elements  connected 
together  at  adjacent  attachment  points  - called  nodes  or  grid  points. 
For  the  case  at  hand, the  elements  will  be  considered  as  short  beams 
connected  end  to  end  to  form  the  longer  beam  structure  of  interest. 

One  then  seeks  to  relate  (at  least  approximately)  the  displacement  at 
any  point  interior  to  an  element  solely  in  terms  of  certain  generalized 
displacements  assumed  at  its  attachments. 

Figure  la  shows  the  beam  structure  of  interest  (corresponding  to 
equation  (8))  broken  down  into  three  shortor  beam  segments  or  eler  ts. 
The  generalized  displacements  at  the  points  of  attachment  consist  of 
one  translation  and  one  rotation.  (We  could  define  more).  It  is  ob- 
vious that  when  adjacent  elements  are  connected  the  element  displace- 
ments at  each  point  of  attachment  must  agree,  i.e.,  must  be  continuous. 

1 2 

e.g.  u3  « u2  , etc. 

This  continuity  requirement  between  element  displacements  therefore 
reduces  the  number  of  generalized  grid  jsoint^  For  the 

three  beam  elements  shown  in  fig  1(a),  the  number  of  independent  dis- 
placements is  thus  reduced  from  twelve  to  eight  upon  connection  of 
the  elements  as  shown  in  figure  1(b). 
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Letting  y(x,t)  represent  the  transverse  displacement  of  the  beam  con- 
tinuum between  the  end  points  of  a given  element  (fig  1(c))  the  pro- 
cedure is  then  to  approximate  y(x,t)  as  a linear  function  of  the 

0 

generalized  element  displacements  u^.  (The  assumption  of  linearity 

e 

will  lead  to  linear  differential  equations  in  the  u^)  Since  there 
e 

are  four  u^  per  element  - all  of  which  are  as  yet  arbitrary  - we  can 
try  a polynomial  expression  with  four  arbitrary  constants , i.e,  a 
cubic: 

y(x)  ■ aQ  ♦ ajX  + a2x  ♦ a3x  ■ ajX 
where,  in  the  latter  notation  we  intend  that  a sum  be  performed  over 
repeated  index  i,  with  i ranging  from  0 to  3. 


The  a^  can  easily  be  determined  from  th*  four  conditions: 
y(o,t)  ■ u*  ; y' (o,t)  ■ u®  ; 

y(t,t)  - u®  ; y’(*,t)  ■ u*  ; 

The  result  can  be  expressed  as  the  vector  product: 

y(x,t)  ■ a(x)u®  ; 0 < l ...  (9) 

explicitly, 

y(x,t)  - { i -3C2*253|  t(C-2C2H3)|  3C2-253|U3-C2)t  > 

» 

where  5 ■ x/l 
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Equation  (9)  constitutes  a formal  discretization  of  the  continuum  in 
that  all  interior  displacement  information  has  been  referred  to  the 
end  point,  or  element,  displacements.  The  goal  then  is  to  determine 
these  displacements  u®(t). 

When  subjected  to  sets  of  applied  forces  and  constraints  (boundary 
conditions),  the  beam  element  responds  according  to  the  laws  of 
mechanics  from  which  we  are  free  to  choose  any  one  of  several  in 
which  to  formulate  the  equations  of  motion  of  the  element.  It  is  very 


convenient  to  employ  the  principle  of  virtual  work  for  dynamic  load- 
ing 5 - which  states  that  in  a virtual  displacement  6y(x)  , of  the 
beam  element  from  its  instantaneous  state  of  equilibrium,  the  incre- 
ment in  strain  energy,  i.e,,  the  virtual  strain  energy,  is  equal  to 
the  sum  of  the  virtual  work  done  by  all  the  forces  including  the 
inertia  loads. 


6U  ■ 6W  - $6yyAdx 


...(10) 


where  6U  represents  the  virtual  elastic  strain  energy  resulting  from 
the  virtual  work  6W  of  the  applied  forces  and  the  virtual  work  of 
the  inertia  forces,  (p  is  the  material  density  and  A the  beam  cross- 
sectional  area.)  In  general  the  elastic  strain  energy  due  to  a vir- 
tual displacement  can  be  written: 


6U  ■ /o5e dxdydz 

where  o is  the  induced  stress  due  to  virtual  strain  6c. 


...(11; 


(5)  Przmieniecki,J,S, , Theory  of  Matrix  Structural  Analysis.  McGraw 
Hill,  1968,  p.  267 


• • ♦ § jg__j|  + n 


The  virtual  work  due  to  applied  forces  f(x,t)  per  unit  length: 


<$W  * /l(x,t)6ydx 
0 

Thus  (10)  becomes: 

A a.. 

/ofiedxdydz  « /f(x,t)6ydx  -ZpySyAdx  ...(12) 

v 0 0 


Hooke's  Law  specifies: 
a ■ 0c 

In  beam  theory,  S «■  E,  Young's  modulus,  and  c * hy",  h being  the 
distance  from  the  beam  neutral  axis  to  the  fiber  in  which  a is  being 
defined.  Substituting  in  (12),  the  left  term  becomes: 

& 4 
/ / Eh2y"«y"dxdydz  - El/  y"6y"dx 
0 A 0 

so  that  (12)  becomes: 

l l 4 .. 

EI/y"6y’'dx  ■ /f(x,t)6ydx  -/  oySyAdx  ...(13) 

0 0 0 

Making  use  of  the  approximation  (9),  i.e,,  y «*  £(x)u® ; 

6y  » £(x)6ue  and  6yM  » a"(x)6ue 
Substituting  these  expressions  in  (13): 

4 4.  . 4_e 
EI/6uea"a”u8dx  * /djjeaf(x,t)dx  - oA/6u  3auedx 


Since  the  virtual  displacements  u®  are  arbitrary: 


A £ 

{ pA/aadx  }iae  ♦ { Elf  aJV'dx  Ju^ 


/af(x,t)dx 

Q— 


...(14) 


(a  bar  over  a quantity  denotes  its  transpose) 


The  nxn  matrix  coefficient  expressions  (n  » 4 for  the  problem  at  hand) 
of  ii  and  u deserve  to  be  called  m®  and  £e  respectively  and  the 
right  hand  term  is  the  force  vector  f*  whose  elements  replace  the 
distributed  and  applied  forces  present  in  the  continuum  problem. 

These  forces  are  considered  as  being  applied  to  the  ends  of  the 
element  e,  The  matrices  m®  and  have  been  evaluated  many  times 
in  the  literature6  and  will  simply  be  repeated  here  for  beam  elements 
of  four  degrees  of  freedom. 


k® 


m 


EI/A3  * 


pAA/420 


12 


61 

412 

symm 

-12 

-61 

12 

61 

21* 

•61 

4.2 

— 

— 

1S6 

221 

4t2 

symm 

54 

131 

156 

-111 

-3t2 

-221 

4t2 

(6)  Przmieniecki , J.S.,  Theory  of  Matrix  Structural  Analysis,  McGraw 
Hill,  1968,  p.  81,  297 
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To  form  the  force  vector  f®,  the  expression  f(x,t)  from  equation  (7) 
is  substituted  into  the  right  side  of  equation  (14), 

l 

thus  f = -fa  m (y  * 2vy'  ♦ g + v yu)$(x-vt)dx 

— 0~  P 

0 

substituting  the  relation  y(x)  = au  : 

i , 

f6  = -m  / a(x)  {aue  ♦ 2v^'u>®  +g  + vV'ue}6(x-vt)dx 

» -mpfcjCt)!!®  + c2(t)u*  * CjftJvi®}  -nipga(vt) 

,,,  0 <vt<£ 

■ 0 otherwise 


where 


Ci(t)  ■ ^(vt)a^(vt) 

c2(t)  ■ 2v£(vt)£' (vt) 

Cj(t)  « v2£(vt)aV(vt) 


Equation  (14)  can  then  be  written: 

(m®  ♦ mpCl(t))ue  + * (i®  * Wpc3 (*))«*  * 

, , , 0<_  vt <1 


or  simply, 


ue(t)ue  ♦ ^(t)u_®  ♦ <,(t)ii®  ■ £®  (t)  ,.,(15) 


14 
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Outside  the  interval  0 <vt«,  , the  time  variant  coefficients  c^(t)  and 
«i(vt)  must  be  replaced  by  zeros.  This  can  be  accomplished  by  nullify- 
ing the  m factor  outside  this  interval. 

P 

Equation  (15)  represents  a set  of  n ordinary  differential  equations 
with  time  dependent  coefficients.  They  are  the  differential  equations 
of  motion  for  any  beam  element  of  density  p and  section  modulus  El 
as  well  as  cross-sectional  area  A.  All  of  these  element  properties 
may  differ  from  element  to  element, 

STRUCTURE  EQUATIONS  OF  MOTION 

The  equations  of  motion  for  the  combined  structure,  i.e,, 

•• 

M(t)U  ♦ C(t)U  ♦ £(t)U  - F(t)  ...(16) 

are  formed  as  follows. 

Each  term  of  equation  (15)  constitutes  a force  - , at  the  element 

attachment  points,  i ■ i thru  n.  When  all  elements  are  joined  the 
resultant  force  at  the  connections  (grid  points)  is  the  sum  of  the 
individual  forces  at  the  attachments.  For  example,  the  inertia  forces 
at  the  right  end  of  the  first  element  (cf,  fig  1)  are  to  be  added  to 
those  at  the  left  end  of  the  second  element  to  yield  the  total  inertia 
force  at  the  connection  grid  point, 

i.e,,  the  inertia  forces  acting  on  the  first  element  - element  #1 
in  figure  1 - in  the  Uj,  u2,  u3,  and  u4  directions  are  the  vector  com- 

1 "1 


ponents : 


1 


Summation  over  repeated  subscripts  is  intended  and  i § j range  from  1 
to  n.  Similarly  for  element  #2: 


<>?  * y?.u .? 
A j 


1 2 

Upon  joining  these  two  elements,  ♦ combines  with  the  force  and 


with  i $ so  that  the  resulting  forces  on  the  two-beam  substructure 


are,  for  the  case  n » 4: 

r .i  n r 


^ * *2 
3 1 


,1  2 

4 ” e 


4 2 

.,.2 


1 

, 1 

1 

, 1 

y 

U 

y 

11 

12 

13 

14 

1 

1 

1 

1 

y 

y 

y 

21 

22 

23 

24 

1 

1 

1 

1 

31 

U33 

U34 

1 

1 

1 1 

1 
« i 

1 

41 

42 

43 

W44 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

11 

2 

U12 

2 

W13 

2 

U14 

2 

. 

21 

2 

W22 

2 

P23 

2 

U24 

2 

31 

y2 

32 

y2 

33 

y2 

34 

2 

, 

41 

2 

W42 

2 

W43 

2 

U44 

Enforcing  the  equality  of  the  displacements,  velocities,  and  acceler* 
ations  of  the  attachments,  i.e. 

••1  "2  ::  . "1  "2 


where  the  upper  case  letter  denotes  that  Uj  , are  grid  point 
accelerations  in  conformance  with  figure  1(b). 


tv 

x 

. a*  ‘j 


• ill 


• • • 


Hence  in  terms  of  grid  point  notation: 


U12 

/ u1 

13  W14 

0 

M 

(. 

M1  U1 

23  W24 

0 

0 

(W33+Wll)tW34+U12) 

u2 

13 

*S« 

* 

Cli43+U21J  (y44+U22J 

u2 

m23 

A, 

0 

2 2 
W31  *32 

*33 

•5. 

0 

u2  u2 

41  *42 

W43 

4 

= M(t)£(t)  the  inertia  forces  acting  on  the  grid  points 

of  a 2-element  substructure.  The  other  forces  in  the  equations  of  ao» 
tion  for  this  structure  are  formed  by  similar  superposition.  Thus  the 
structure  equations  of  motion  are  formed  by  ***• 

element  matrices  wherever  a grid  point  connection  is  made.  For  an  N • 
element  beam  there  will  be  N-l  such  overlaps  (shown  schematically  in 
figure  2).  Each  overlap  will  contain  n/2  entries  from  the  lower  right 
corner  of  the  nxn  matrix  corresponding  to  the  element  to  the  left  of 
the  grid  point  - and  n/2  entries  from  the  upper  left  comer  of  the 
matrix  corresponding  to  the  element  to  the  right  of  this  grid  point. 
These  overlapping  elements  are  to  be  summed. 


• • • • • • 
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Equations  of  Motion  for  the  Moving  Mass  Problem 


It  is  to  be  noticed  from  equation  (15)  that  in  general,  the  coeffi» 
cient  matrices  v£  , yf , consist  of  a constant  part  (null  in  the 

case  of  ^®)  and  a time  variant  part  which  derives  from  the  moving  mass 
nip.  These  time  variant  elements  are  nul?  except  when  t is  such  that  nip 
is  located  within  the  length  of  a particular  beam  element.  Only  then 
are  the  time-variant  components  of  the  corresponding  element  matrix 
finite.  Thus  figure  2, one  conceives  of  a conventional  matrix  of 
constant  coefficient  multipliers  of  the  acceleration,  velocity,  and 
displacement  terms  plus  a time  variant  set  of  components  which  propagate 
in  a band  along  the  diagonal  of  each  structure  coefficient  matrix  as  the 
moving  mass  traverses  the  beam  in  time.  Thus  at  any  instant  only  n of 
the  structure  equations  of  motion  possess  time  variant  coefficients  - 
n being  the  number  of  element  displacements  (degrees  of  freedom)  con- 
sidered for  each  beam  element  • four  for  the  case  at  hand.  Thus  the 
prospect  of  solving  the  full  set  of  equations  numerically  - without 
incurring  extraordinarily  long  computer  run  times  - would  appear  to 
be  good.  For  example  it  is  not  uncommon  to  solve  via  computer,  a 
fifty  degree  of  freedom  transient  problem  in  conventional  structure 
dynamics  via  the  finite  element  technique  - whert  all  of  the  coefficients 
are  constant  in  time.  It  should  therefore  involve  only  a moderate  in- 
crease in  computation  time  to  allow  four  of  these  equations  to  take 
on  time- variant  coefficients  as  in  our  problem  of  the  moving  concen- 
trated mass.  Roughly  speaking  one  might  expect  that  each  time-variant 
matrix  element  will  create  additional  computation  no  greater  than  that 
caused  by  adding  another  degree  of  freedom  to  a conventional  constant 
coefficient  problem.  Thus  a 50  degree  of  freedom  problem  - e.g.. 


twenty- four  connected  beam  elements  subjected  to  a moving  concentrated 
mass  - could  be  solved  with  a computation  time  not  in  excess  of  a 
ninety-eight  degree  of  freedom  problem  in  which  all  the  matrix  coef- 
ficients are  constant.  ( yfCt),  £?t),  and  £e(t)  each  comprise  sixteen 
time  dependent  components). 

Boundary  Conditions 

A great  convenience  of  the  finite  element  procedure  - as  compared 
say,  to  the  Ritz  or  Galerkin  procedures  - is  that  all  ambiguity  is 
removed  in  choosing  the  boundary  conditions  to  enforce.  ( This  is  due 
to  the  particular  stage  of  deduction  at  which  the  finite  element  ideal- 
ization is  invoked  in  a variational  procedure*).  In  practice  all  one 
has  to  do  is  mimic  physical  reality.  For  example,  a beam  with  hinge 
(simple)  supports  at  each  end  requires  that  the  corresponding  displace- 
ments vanish;  e„g.,  in  a three-element  beam  model,  Uj  » U7  ■ 0.  Simi- 
larly a clamped  cantilevered  beam  would  insist  that  Uj  ■ U2  ■ 0,  In- 
stead of  specification  of  particular  zero  values,  however,  it  is  more 
efficient  to  merely  delete  the  corresponding  rows  and  columns  from  the 
coefficient  matrices  M , £,  and  £ of  equation  (16).  Thus  for  the  case 
of  a three-element  hinged-hij  d beam,  we  simply  delete  the  first  and 
the  seventh  rows  and  columns  •-  om  these  matrices.  Similarly  we  delete 
the  corresponding  elements  from  any  force  vector  £ , appearing  on  the 
right  hand  side  of  this  equation. 


‘Conversations  with  Dr.  Gary  Anderson,  Applied  Mathematics  and  Mechanics 
Div,  Benet  Weapons  Laboratory,  Watervliet  Arsenal,  Watervliet,  N.Y. 
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PROBLEM  STATEMENT  AND  SOLUTION 

For  the  purpose  of  comparing  results  achieved  by  finite  elements 
with  those  appearing  in  the  literature,  the  problem  to  be  solved  is  one 
of  three  beam  elements  connected  end  to  end  to  form  a simply  supported 
beam  having  uniform  cross  section  and  material  properties.  A moving 
mass,  nip  , is  assumed  to  traverse  the  beam  from  left  to  right  at  a 
constant  velocity  v - see  figure  3, 

Deleting  the  first  and  seventh  rows  and  columns  from  the  coefficient 
matrices  N[,£  and  £ along  with  the  first  and  seventh  components  of  the 
vector  £ of  equation  (16)  results  in  a six  degree  of  freedom  problem 
- that  is,  six  coupled  equations  of  forced  motion.  These  equations 
will  be  solved  numerically  with  nip  and  its  velocity  v,  serving  as 
parameters . 

In  reality  we  have  no  clairvoyance  to  guide  the  choice  of  values 
for  nip  and  v except,  of  course,  to  repeat  those  used  in  the  literature 
so  that  a comparison  may  be  made.  It  appears,  however,  that  the 
values  chosen  by  Ayre,  et  al,  were  not  completely  arbitrary  in  that 
certain  values  of  v will  cause  resonant  (secular)  behavior  in  the 
moving  force  problem.  From  equation  (3)  we  can  verify  that  resonance 
will  indeed  occur  for  any  value  v ■ v*,  such  that  m^(v*)  ■ 0.  Al- 
though applicable  only  to  moving  force  problems,  it  is  intuitively 
plausible  that  extraordinary  behavior  in  the  moving  mass  problem 
might  occur  for  values  of  velocity  not  far  removed  from  these  values 
v*.  Having  no  'closed  form'  analytical  solution  with  which  to  an- 
ticipate points  of  singular  • or  otherwise  interesting  - behavior,  (me 
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can  do  no  better  than  to  at  least  allow  the  velocity  of  the  mass  to 
range  through  the  lowest  v*  which  will  cause  m^  to  vanish.  This 
apparently  was  the  reasoning  behind  the  choice  of  velocities  investiga- 
ted by  Ayre,  et  al,  who  performed  experiments  in  which  v was  chosen 
as  v*/4,  v*/2,  and  v*  for  the  moving  force  and  v*/4  and  v*/2  for  the 
moving  mass  - v*  being  the  first  value  of  v to  cause  m^  to  vanish, 
i.e.,  the  first  resonant  or  'critical'  velocity  of  the  moving  force 
solution.  (Maintaining  contact  between  the  beam  and  mass  for  velocities 
higher  than  v*/2  was  apparently  impossible.) 

For  the  simply  supported  beam  being  considered: 

mii(v)  * aiv2  " “i  ...(17) 

where  - i*/L,  L being  the  overall  length  of  the  contained  beam 

ft 

structure  and  ■ a*  /EI/oA 

Table  1 lists  t)  e material  constants  and  dimensions  employed  for  the 
three  identical  beam  elements  used  in  this  work.  From  these  values 
one  calculates  from  (17)  for  i«l: 

v*  ■ 899,13  in/sec 

We  intend  to  examine  numerical  solutions  for  the  moving  force  and 

* 

the  moving  mass  throughout  the  range  0 <v  ^v  , 

Two  primary  references  will  be  used  as  basis'  of  comparison: 

O 

(i)  The  exact  solution  for  the  moving  force  solution  which  for 


(8)  Nowacki,  W,  Dynamics  of  Elastic  Systems,  Chapman  5 Hall  Ltd, 
London,  1963,  p,  lib 
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TABLE  1.  MATERIAL  PROPERTIES  AND  NOMENCLATURE 


pAL 


Materi al 


Density 


Overall  Beam  Length 


Length  of  each  element 


Beam  Mass 


aluminum 


3,14  x 10"4lbsec2/in4 


3Z 

120,0  in 


m Moving  Mass 

P 


E Young's  Modulus 


Beam  Cross  Sect.  Aroa 


0.0,  pAL/4,  pAL/2 
1.0  x 107  psi 
31 


I Beam  Area  Moment  of  Inertia  (fixed  by  choice  of  wl  below) 


h . Beam  Thickness 
1^/2*  Beam  Fundamental  Frequency 

v Mass  Velocity 


2.0  in. 


1.25  hz 


100.,  v*/4,  v*/2,  v* 


Fundamental  Resonant  Velocity  899,13  in/sec 


• • • • 
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the  case  of  a moving  concentrated  downward  force  -mpS(x-vt)  is  : 

y(x,t)  » -2mpg/oAL  Jsin^xCcinVSinu^t  - «nsinanvt)/un(a2v2-o>J) 

... CiS) 

The  finite  element  solutions  for  the  moving  mass  problem  (Equation  (16) 
after  imposing  support  constraints)  with  • 0 will  be  compared  with 
computations  of  (18)  above, 

2 

(ii)  The  experimental  work  of  Ayre,  Jacobsen  and  lisu  , The  mass 

• * 

velocities  employed  in  this  work  were  quasi-static  v ■ 0,  and  v ■ v /4, 
v*/2,  and  v*.  The  moving  mass  values  chosen  were  nip  ■ 0,  and  nip  ■ 
pAL/4  and  pAL/2  where  pAL  is  the  total  beam  mass. 


Equation  (16)  is  solved  using  Hamming's  modified  predictor-corrector 
method  which  uses  fourth  order  Runge-Kutta  method  suggested  by  Ralston10 
for  adjustment  of  the  initial  increment  and  for  computation  of  starting 
values r The  method  is  taken  directly  from  the  IBM  Scientific  Subroutine 
Package  for  the  IBM  System  #360,  re:  Programmer's  Manual  # H20-0205 
available  this  laboratory.  The  method  was  found  to  be  about  four  times 
faster  than  using  the  Runge-Kutta  method  throughout  the  entire  problem. 
In  general,  run  times  in  the  order  of  20  - 30  minutes  on  the  IBM  model 


44  computer. 


(2)  Ayre,  R.S.,  Jaccbsen,  L.S.,  and  Hsu,  C.S.,  "Transverse  Vibration 
of  One  and  of  Two  Span  Beams  Under  the  Action  of  a Moving  Mass  Load," 
Proc.  of  First  National  Congress  on  Applied  Mechanics,  June  1951 


(10)  Ralston  and  Wilf,  Mathematical  Methods  for  Digital  Computers, 
Wiley  and  Sons,  New  York,  London,  1960,  pp.  95-109 
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Results  - Moving  Force 


Figures  4 through  7 show  the  transverse  displacement  of  the  beam  at  the 


grid  point  locations  x » L/3,  2L/3  as  computed  either  by  the  finite 


element  method  or  by  evaluation  of  the  exact  solution  (18)  - any  dis- 


crepancy between  the  results  being  too  small  to  be  discerned  even  in 


plots  of  this  scale.  This  is  substantially  less  error  than  appears  in 

4,9 

other  treatments  yielding  approximate  results. 


Using  the  relationship  (9),  i.e,,  y(x,t)  ■ £(x)i je(t)  with  x ■ vt, 


gives  the  displacement  directly  beneath  the  moving  load  to  be  compared 


with  the  results  of  Ayre,  et  al,  who  recorded  displacement  information 


exclusively  at  this  location.  Figures  8 • a,b,c  shew  these  comparisons 


and  it  is  obvious  that  the  agreement  with  experiment  is  much  better  in 


figures  b and  c than  in  a.  Actually,  Ayre  and  his  co-workers  experienced 


considerable  experimental  difficulty  when  the  force  was  translated  at 


v /4,  Quoting  from  their  publication  in  which  the  authors  remark  on 


their  disagreement  with  Schallenkamp's  theoretical  solution  for  the 


moving  force  (evidently  a three  term  approximation  of  expression  (18)): 


..."The  agreement  is  generally  good  [except  at  v /4]  where  it  has  been 
found  that  comparatively  small  errors  in  velocity  may  result  in  marked 


differences  in  the  shape  of  the  trajectory,"  ...  In  8a  the  theoreti- 


cal curve  used  by  Ayre  as  a basis  of  comparison  has  been  included. 


(2)  Ayre,  R.S.,  Jacobsen,  L.S,,  and  Hsu,  C.S.,  "Transverse  Vibration 
of  One  and  of  Two  Span  Beams  Under  the  Action  of  a Moving  Mass  Load," 
Proc.  of  First  National  Congress  on  Applied  Mechanics,  June  1951 


(4)  Schallenkamp,  A.,  "Schwingungen  von  Tragem  bei  bewegten  Lasten," 
Ingenieur-Archiv,  v.8,  1937,  pp.  182-198 


(9)  Hutton,  D.V.,  and  Counts,  J,,  "Deflections  of  a Beam  Carrying  a 
Moving  Mass,"  Trans.  ASME,  Sept,  1974,  p,  803 
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Figure  6 
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Figure  7 
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Transverse  Motijn  at  Grid  Points  in  Response  to  Moving  Force 
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The  exceptional  agreement  (using  only  three  finite  beam  elements)  of 
the  finite  element  results  for  the  moving  force  gives  confidence  in  ex- 
tending the  technique  to  the  moving  mass  problem  - for  which  no  simple 
solution  as  (18)  exists. 

Results  - Moving  Mass 

Figures  9 through  14  show  the  transverse  displacement  of  the  beam  at  the 
grid  point  locations  x » L/3,  2L/3  as  computed  by  the  finite  element 
procedure  for  various  cases  of  velocity  and  mass  values.  Again  use 
will  be  made  of  relationship  (9)  to  convert  this  information  to  dis- 
placement beneath  the  moving  mass  so  that  a direct  comparison  with  the 
experimental  results  of  Ayre,  Jacobsen,  and  ilsu  can  be  made.  Figure 
15  - a,b,c  shows  this  comparison  to  be  quite  good  - excellent  agree- 
ment occurring  at  the  grid  point  locations,  A closer  look  shows,  however, 
that  the  slopes  of  the  curves  generated  by  the  finite  element  analysis 
are  discontinuous  at  the  grid  point  locations.  The  reason  for  this  is 
that  the  displacement  approximation  (9)  is  built  from  cubic  polynomials 
wnich  are  not  continuous  in  the  second  derivative  at  the  grid  point 
connections.  Indeed,  only  continuity  in  y and  y'  were  demanded  in 
constructing  these  p .omials. 

The  discontinuities  of  y" (x)  at  grid  points  might  not  be  serious  if 

it  was  not  for  a 'force'  in  the  continuum  equation  of  motion  (8),,, 

mpv  y"(x).  Thus  at  higher  velocities  jump  discontinuities  in  y"  will 

cause  increasingly  powerful  disturbances  which  are  nonphysical  in 

character.  This  is  apt  to  be  especially  influential  in  armament 

"Note:  All  displacements  are  normalized  with  respect  to  the  displace- 
ment which  occurs  at  midspan  due  to  a static  load  at  this  point. 
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Figure  12 
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Figure  15.  Transverse  Motion  of  Moving  Mass 


applications  involving  velocities  considerably  larger  than  v*. 

One  approach  toward  rectification  of  this  problem  requires  tKe 
formation  of  higher-ordered  polynomials  which  are  continuous  in  the 
second  derivative  across  the  grid  points  of  the  structure.  Higher 
polynomials  have  additional  coefficients  and,  as  in  the  previous  deri- 
vation, these  coefficients  are  equal  in  number  to  the  element  degrees 
of  freedom,  i.e.,  the  generalized  displacements  u®  , 

The  establishment  of  polynomial  displacement  functions  proceeds  in 
the  same  manner  as  before,  with  the  added  conditions: 

y"(0)  ■ u®  ; y"(fc)  ■ uj 

The  result  is  that  y(x)  will  be  approximated  interior  to  an  element  by: 

y(x)  « a(x)u®  ; 

where  £ and  u%re  now  vectors  with  six  components  instead  of  four. 

The  new  (*)  element  matrices  are  calculated  according  to  equation 
(14).  The  multiplication  and  integration  was  accomplished  analytically 
by  a computer  program  called  MANIP  which  was  written  in  the  FORMAC 
language  • giving  analytical  expressions  as  output.  Similarly,  new 
expressions  for  c^,  c2,  and  c^  for  use  in  equation  (IS)  were  formed 
using  this  program,  and  m compose  the  constant  elements  of  the 
matrices  and  j^e  and  are  found  to  be: 
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Figure  16  shows  the  improved  results  obtained  through  the  use  of  the 
quintic  polynomial.  Comparison  should  be  made  with  the  results  obtained 
via  the  cubic  polynomial  as  shown  in  figure  15(a).  The  price  paid, 
however,  is  increased  computation  time  due  to  the  additional  degrees 
of  freedom  induced  through  the  use  of  the  quintic.  One  also  faces 
increasing  amount  of  ill-conditioning  in  the  matrix  equations  of  motion 
due  to  the  widely  differing  magnitudes  induced  in  the  coefficient  ma- 
trices and  their  corresponding  output  variables.  Future  work  will 
concentrate  on  these  problems  and  the  use  of  less  time-consuming 
integration  algorithms. 
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Figure  16.  Improved  Results  Using  Quintic  Shape  Function. 
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